	/*This program computes logit demand and simulates corresponding price changes assuming Bertrand competition with product differentiation*/

	#delimit;
	set matsize 10000;
	clear;
	clear matrix;
	mata: mata clear;
	local nbrands=9;
	set seed 070910432;

	*quietly{;
		set mem 250m;
		set more off;
		do "c:\data\restat_12_9\programs\csolve_do\csolve";
		do "c:\data\restat_12_9\programs\csolve_do\callf";
		do "c:\data\restat_12_9\programs\csolve_do\callf_setup";
		do "c:\data\restat_12_9\programs\csolve_do\\antilogit_7_8";
		/*get top level elasticity*/
		mata: mata matuse "c:/data/restat_12_9/analysis_data/post_top_level";
		
		mata: st_matrix("ols_top_b",ols_b_top[1,1]);
		mata: st_matrix("iv_top_b",iv_b_top[1,1]);
		
		mata: st_matrix("ols_top_mean",ols_b_top);
		mata: st_matrix("iv_top_mean",iv_b_top);
		mata: st_matrix("ols_top_vc",ols_vc_top);
		mata: st_matrix("iv_top_vc",iv_vc_top);

		
		
		scalar ols_top_b=ols_top_b[1,1];
		mata: ols_top_b=ols_b_top[1,1];
		
		scalar iv_top_b=iv_top_b[1,1];
		mata: iv_top_b=iv_b_top[1,1];
		
		*scalar top_b=-1.15;
		
		
		use "c:\data\restat_12_9\analysis_data\post_hl_syruplogitdata";
		qui tab BRAND2, gen(bra);
		qui tab REGION, gen(reg);
		qui tab MONTH, gen(m);
		qui tab t, gen(t);
		local reglist "reg2";
		forval x=3/47{;
			local reglist "`reglist' reg`x'";
		};
		
		su AVOL;
		scalar vt=r(N)*r(mean);
		forval x=1/`nbrands'{;
			su AVOL if bra`x'==1;
			scalar vshare`x'=r(N)*r(mean)/vt;
		};

		forval x=1/`nbrands'{;
			su price if bra`x'==1;
			scalar price`x'=r(mean);
		};
		scalar wavgp=vshare1*price1+vshare2*price2+vshare3*price3+vshare4*price4+vshare5*price5+vshare6*price6+vshare7*price7+vshare8*price8+vshare9*price9;
		
		forval x=1/`nbrands'{;
			forval y=2/47{;
				gen b`x'r`y'=reg`y'*bra`x';
			};
		};
		
		forval x=1/`nbrands'{;
			forval y=2/12{;
				gen m`x'r`y'=bra`x'*m`y';
			};
		};
		
		forval x=1/`nbrands'{;
			gen trend`x'=bra`x'*t;
		};
		
		
		forval x=1/8{;
			forval y=2/47{;
				local brandxreg "`brandxreg' b`x'r`y'";
			};
		};
		
		
		forval x=1/8{;
			forval y=2/12{;
				local brandxmonth "`brandxmonth' m`x'r`y'";
			};
		};
		
		
		
		forval x=1/8{;
			local trends "`trends' trend`x'";
		};
		di "`trends'";
		
	
		
	
		local logit "price_plprice bra1 bra2 bra3 bra4 bra5 bra6 bra7 bra8 `reglist' t2 t3 t4 t5 t6 t7 t8 t9 t10 t11 t12 t13 t14 t15 t16 t17 t18 t19 t20 t21 t22 t23 t24 
		t25 t26 t27 t28 t29 t30 t31 t32 t33 t34 
		t35 t36 t37 t38 t39 t40 t41 t42";
		
		local logit2 "price_plprice bra1 bra2 bra3 bra4 bra5 bra6 bra7 bra8 `brandxreg' `brandxmonth' ";
		
		sort BRAND REGION t;
		
		
		sort BRAND REGION t;
		by BRAND REGION: gen T=_n;
		sort T REGION BRAND;
		bysort T: gen id=_n;
		tsset id T;
		
		noisily xi: newey2 lshare_plshare `logit'  if BRAND2~="PRIVATE LABEL", lag(1) nocons;
		noisily xi: newey2 lshare_plshare `logit2'  if BRAND2~="PRIVATE LABEL", lag(1) nocons;
		
		matrix ols_logitB=e(b);
		matrix ols_logitVC=e(V); 
		scalar alpha=_b[price_plprice];
		matrix olsb=_b[price_plprice];
		
		
		forval x=1/`nbrands'{;
			forval y=1/`nbrands'{;
				sca olse`x'b`y'=-price`x'*alpha*vshare`x'+ols_top_b*vshare`x'*price`x'/wavgp;
			};
		};
		
		forval y=1/`nbrands'{;
			sca olse`y'b`y'=price`y'*alpha*(1-vshare`y')+ols_top_b*vshare`y'*price`y'/wavgp;
		};
		
		mat olselas=J(9,9,0);
		forval x=1/`nbrands'{;
			forval y=1/`nbrands'{;
				mat olselas[`x',`y']=olse`x'b`y';
			};
		};
		mat li olselas;
		
		
		noisily xi: newey2 lshare_plshare bra1 bra2 bra3 bra4 bra5 bra6 bra7 bra8 `reglist' t2 t3 t4 t5 t6 t7 t8 t9
		t10 t11 t12 t13 t14 t15 t16 t17
		t18 t19 t20 t21 t22 t23 t24 t25 t26 t27 t28 t29 t30 t31 t32 t33
		t34 t35 t36 t37 t38 t39 t40 t41 t42 (price_plprice=tiprice1-tiprice42) if BRAND2~="PRIVATE LABEL", lag(1) nocons;

		noisily xi: newey2 lshare_plshare bra1 bra2 bra3 bra4 bra5 bra6 bra7 bra8 `brandxreg' `brandxmonth'  (price_plprice=iprice) if BRAND2~="PRIVATE LABEL", lag(1) nocons;
		
		matrix iv_logitB=e(b);
		matrix iv_logitVC=e(V);
		
		scalar ivalpha=_b[price_plprice];
		matrix ivb=_b[price_plprice];
		
		
		forval x=1/`nbrands'{;
			forval y=1/`nbrands'{;
				sca ive`x'b`y'=-ivalpha*vshare`x'*price`x'+iv_top_b*vshare`x'*price`x'/wavgp;
			};
		};
		
		forval y=1/`nbrands'{;
			sca ive`y'b`y'=ivalpha*price`y'*(1-vshare`y')+iv_top_b*vshare`y'*price`y'/wavgp;
		};
		mat ivelas=J(9,9,0);
		forval x=1/`nbrands'{;
			forval y=1/`nbrands'{;
				mat ivelas[`x',`y']=ive`x'b`y';
			};
		};
		mat li ivelas;
		


		mat vsbar=J(9,47,0);
		mat pbar=J(9,47,0);
		forval r=1/47{;
			sum AVOL if reg`r'==1;
			sca vt`r'=r(N)*r(mean);
		};
	
	
	
		forval r=1/47{;
			forval b=1/9{;
				qui su AVOL if reg`r'==1&bra`b'==1;
				mat vsbar[`b',`r']=r(N)*r(mean)/vt`r';
				qui su price if reg`r'==1&bra`b'==1;
				mat pbar[`b',`r']=r(mean);
			};	
		};		

		
		
	#delimit cr;
		mata:
			pbar=st_matrix("pbar")
			vsbar=st_matrix("vsbar")
			
			
			pre_m=I(`nbrands')
			pre_m[1,2]=1
			pre_m[2,1]=1
			pre_m[3,4]=1
			pre_m[4,3]=1
			pre_m[5,6]=1
			pre_m[6,5]=1
			pre_m[7,8]=1
			pre_m[8,7]=1
			pre_m[9,9]=1
			
			
			post_m=pre_m
			post_m[5,7]=1
			post_m[5,8]=1
			post_m[6,7]=1
			post_m[6,8]=1
			
			post_m[7,5]=1
			post_m[7,6]=1
			post_m[8,5]=1
			post_m[8,6]=1
			
			olslogit=J(`nbrands',47,0)
			ivlogit=J(`nbrands',47,0)
			ols_solver_flag=J(47,1,0)
			iv_solver_flag=J(47,1,0)
			olsb=st_matrix("olsb")
			ivb=st_matrix("ivb")
			olsb=-olsb
			ivb=-ivb
			ols_e=-ols_top_b
			iv_e=-iv_top_b
		
			
			for(i=1;i<=cols(pbar);i++){
					pinit=pbar[.,i]
					a=csolve(&anti(),pinit,1e-6,400,vsbar[.,i],pbar[.,i],olsb,ols_e,pre_m,post_m)
					sol=a[1..9,.]
					olslogit[.,i]=sol
					ols_solver_flag[i,1]=a[10,1]
			}
		olslogit=olslogit'
		ans=(olslogit-pbar'):/pbar'
		/*take only markets for which we could compute a solution*/
		ols_sim=100*47/(colsum(J(47,1,1)-!rowsum(ans)))*mean(ans,1)
		
			for(i=1;i<=cols(pbar);i++){
					pinit=pbar[.,i]
					a=csolve(&anti(),pinit,1e-6,400,vsbar[.,i],pbar[.,i],ivb,ols_e,pre_m,post_m)
					sol=a[1..9,.]
					ivlogit[.,i]=sol
					iv_solver_flag[i,1]=a[10,1]
			}
		ivlogit=ivlogit'
		ans=(ivlogit-pbar'):/pbar'
		iv_sim=100*47/(colsum(J(47,1,1)-!rowsum(ans)))*mean(ans,1)
	end		


/***Now that we have our point estimates of costs and approximate simulated price changes, we draw from asymptotic approx of distributions
of estimators and recalculate for each draw.  will look at empirical distribution of these calculated costs and changes to do inference.****/	
#delimit;
drop _all;
/*number of bootstrap iterations*/


local its=1000;
drawnorm `logit2', n(`its') mean(ols_logitB) cov(ols_logitVC);
mkmat _all;

drop _all;
drawnorm e_topbs b2-b60, n(`its') mean(ols_top_mean) cov(ols_top_vc);
mkmat e_topbs;
drop _all;




use "c:\data\restat_12_9\analysis_data\post_hl_syruplogitdata";

tab REGION, gen(reg);
tab BRAND2, gen(bra);


forval x=1/47{;
    mat bpbar`x'=J(`nbrands',`its',0);  
	mat bv`x'=J(`nbrands',`its',0);  
    mat bvsbar`x'=J(`nbrands',`its',0);  
};


forval x=1/47{;
    forval w=1/`its'{;
        preserve;
            keep if reg`x'==1;
            bsample;
            forval y=1/`nbrands'{;
                qui sum price if bra`y'==1;
                matrix bpbar`x'[`y',`w']=_result(3);
                qui sum AVOL if bra`y'==1;
                matrix bv`x'[`y',`w']=_result(3);				
            };
        restore;
    };
};
#delimit cr;
mata:
    bp1=st_matrix("bpbar1")
    bp2=st_matrix("bpbar2")
    bp3=st_matrix("bpbar3")
    bp4=st_matrix("bpbar4")
    bp5=st_matrix("bpbar5")
    bp6=st_matrix("bpbar6")
    bp7=st_matrix("bpbar7")
    bp8=st_matrix("bpbar8")
    bp9=st_matrix("bpbar9")
    bp10=st_matrix("bpbar10")
    bp11=st_matrix("bpbar11")
    bp12=st_matrix("bpbar12")
    bp13=st_matrix("bpbar13")
    bp14=st_matrix("bpbar14")
    bp15=st_matrix("bpbar15")
    bp16=st_matrix("bpbar16")
    bp17=st_matrix("bpbar17")
    bp18=st_matrix("bpbar18")
    bp19=st_matrix("bpbar19")
    bp20=st_matrix("bpbar20")
    bp21=st_matrix("bpbar21")
    bp22=st_matrix("bpbar22")
    bp23=st_matrix("bpbar23")
    bp24=st_matrix("bpbar24")
    bp25=st_matrix("bpbar25")
    bp26=st_matrix("bpbar26")
    bp27=st_matrix("bpbar27")
    bp28=st_matrix("bpbar28")
    bp29=st_matrix("bpbar29")
    bp30=st_matrix("bpbar30")
    bp31=st_matrix("bpbar31")
    bp32=st_matrix("bpbar32")
    bp33=st_matrix("bpbar33")
    bp34=st_matrix("bpbar34")
    bp35=st_matrix("bpbar35")
    bp36=st_matrix("bpbar36")
    bp37=st_matrix("bpbar37")
    bp38=st_matrix("bpbar38")
    bp39=st_matrix("bpbar39")
    bp40=st_matrix("bpbar40")
    bp41=st_matrix("bpbar41")
    bp42=st_matrix("bpbar42")
    bp43=st_matrix("bpbar43")
    bp44=st_matrix("bpbar44")
    bp45=st_matrix("bpbar45")
    bp46=st_matrix("bpbar46")
    bp47=st_matrix("bpbar47")

    ppoint=J(1,47,NULL)
    ppoint[1]=&bp1
    ppoint[2]=&bp2
    ppoint[3]=&bp3
    ppoint[4]=&bp4
    ppoint[5]=&bp5
    ppoint[6]=&bp6
    ppoint[7]=&bp7
    ppoint[8]=&bp8
    ppoint[9]=&bp9
    ppoint[10]=&bp10
    ppoint[11]=&bp11
    ppoint[12]=&bp12
    ppoint[13]=&bp13
    ppoint[14]=&bp14
    ppoint[15]=&bp15
    ppoint[16]=&bp16
    ppoint[17]=&bp17
    ppoint[18]=&bp18
    ppoint[19]=&bp19
    ppoint[20]=&bp20
    ppoint[21]=&bp21
    ppoint[22]=&bp22
    ppoint[23]=&bp23
    ppoint[24]=&bp24
    ppoint[25]=&bp25
    ppoint[26]=&bp26
    ppoint[27]=&bp27
    ppoint[28]=&bp28
    ppoint[29]=&bp29
    ppoint[30]=&bp30
    ppoint[31]=&bp31
    ppoint[32]=&bp32
    ppoint[33]=&bp33
    ppoint[34]=&bp34
    ppoint[35]=&bp35
    ppoint[36]=&bp36
    ppoint[37]=&bp37
    ppoint[38]=&bp38
    ppoint[39]=&bp39
    ppoint[40]=&bp40
    ppoint[41]=&bp41
    ppoint[42]=&bp42
    ppoint[43]=&bp43
    ppoint[44]=&bp44
    ppoint[45]=&bp45
    ppoint[46]=&bp46
    ppoint[47]=&bp47    

    bv1=st_matrix("bv1")
    bv2=st_matrix("bv2")
    bv3=st_matrix("bv3")
    bv4=st_matrix("bv4")
    bv5=st_matrix("bv5")
    bv6=st_matrix("bv6")
    bv7=st_matrix("bv7")
    bv8=st_matrix("bv8")
    bv9=st_matrix("bv9")
    bv10=st_matrix("bv10")
    bv11=st_matrix("bv11")
    bv12=st_matrix("bv12")
    bv13=st_matrix("bv13")
    bv14=st_matrix("bv14")
    bv15=st_matrix("bv15")
    bv16=st_matrix("bv16")
    bv17=st_matrix("bv17")
    bv18=st_matrix("bv18")
    bv19=st_matrix("bv19")
    bv20=st_matrix("bv20")
    bv21=st_matrix("bv21")
    bv22=st_matrix("bv22")
    bv23=st_matrix("bv23")
    bv24=st_matrix("bv24")
    bv25=st_matrix("bv25")
    bv26=st_matrix("bv26")
    bv27=st_matrix("bv27")
    bv28=st_matrix("bv28")
    bv29=st_matrix("bv29")
    bv30=st_matrix("bv30")
    bv31=st_matrix("bv31")
    bv32=st_matrix("bv32")
    bv33=st_matrix("bv33")
    bv34=st_matrix("bv34")
    bv35=st_matrix("bv35")
    bv36=st_matrix("bv36")
    bv37=st_matrix("bv37")
    bv38=st_matrix("bv38")
    bv39=st_matrix("bv39")
    bv40=st_matrix("bv40")
    bv41=st_matrix("bv41")
    bv42=st_matrix("bv42")
    bv43=st_matrix("bv43")
    bv44=st_matrix("bv44")
    bv45=st_matrix("bv45")
    bv46=st_matrix("bv46")
    bv47=st_matrix("bv47")
    vpoint=J(1,47,NULL)
    vpoint[1]=&bv1
    vpoint[2]=&bv2
    vpoint[3]=&bv3
    vpoint[4]=&bv4
    vpoint[5]=&bv5
    vpoint[6]=&bv6
    vpoint[7]=&bv7
    vpoint[8]=&bv8
    vpoint[9]=&bv9
    vpoint[10]=&bv10
    vpoint[11]=&bv11
    vpoint[12]=&bv12
    vpoint[13]=&bv13
    vpoint[14]=&bv14
    vpoint[15]=&bv15
    vpoint[16]=&bv16
    vpoint[17]=&bv17
    vpoint[18]=&bv18
    vpoint[19]=&bv19
    vpoint[20]=&bv20
    vpoint[21]=&bv21
    vpoint[22]=&bv22
    vpoint[23]=&bv23
    vpoint[24]=&bv24
    vpoint[25]=&bv25
    vpoint[26]=&bv26
    vpoint[27]=&bv27
    vpoint[28]=&bv28
    vpoint[29]=&bv29
    vpoint[30]=&bv30
    vpoint[31]=&bv31
    vpoint[32]=&bv32
    vpoint[33]=&bv33
    vpoint[34]=&bv34
    vpoint[35]=&bv35
    vpoint[36]=&bv36
    vpoint[37]=&bv37
    vpoint[38]=&bv38
    vpoint[39]=&bv39
    vpoint[40]=&bv40
    vpoint[41]=&bv41
    vpoint[42]=&bv42
    vpoint[43]=&bv43
    vpoint[44]=&bv44
    vpoint[45]=&bv45
    vpoint[46]=&bv46
    vpoint[47]=&bv47

    b=st_matrix("price_plprice")
	e=st_matrix("e_topbs")

end

mata    
    bslogit=J(1000,`nbrands',0)
	for(z=1;z<=1000;z++){
		bz=-b[z,1]
		ez=-e[z,1]		
        temp=J(47,9,0)

		for(j=1;j<=47;j++){
            	p=(*ppoint[j])[.,z]
	            pinit=(*ppoint[j])[.,z]
				v=(*vpoint[j])[.,z]
				vsbar=v:/colsum(v)
	            a=csolve(&anti(),pinit,1e-6,400,vsbar,p,bz,ez,pre_m,post_m)
      	        
				a=a[1..9,.]
            	a=(a-p):/p
	            temp[j,.]=a'
	      }
	bslogit[z,.]=100*47/(colsum(J(47,1,1)-!rowsum(temp)))*mean(temp,1)
	}
end
mata: st_matrix("bslogit", bslogit)
drop _all
svmat bslogit
save "c:\data\restat_12_9\bootstrap_results\post_ols_syruplogit", replace





#delimit;
drop _all;
/*number of bootstrap iterations*/
local its=1000;
	local reglist "reg2";
	forval x=3/47{;
		local reglist "`reglist' reg`x'";
	};

local logit "price_plprice bra1 bra2 bra3 bra4 bra5 bra6 bra7 bra8 `brandxreg' `brandxmonth'";
drawnorm `logit', n(`its') mean(iv_logitB) cov(iv_logitVC);
mkmat _all;
drop _all;
drawnorm e_topbs b2-b60, n(`its') mean(ols_top_mean) cov(ols_top_vc);
mkmat e_topbs;
drop _all;




#delimit cr;
mata:
    b=st_matrix("price_plprice")
	e=st_matrix("e_topbs")
    bslogit=J(1000,9,0)
	

	for(z=1;z<=1000;z++){
		bz=-b[z,1]
		ez=-e[z,1]		
        temp=J(47,9,0)

		for(j=1;j<=47;j++){
            	p=(*ppoint[j])[.,z]
	            pinit=(*ppoint[j])[.,z]
	            v=(*vpoint[j])[.,z]
				vsbar=v:/colsum(v)
	            a=csolve(&anti(),pinit,1e-6,400,vsbar,p,bz,ez,pre_m,post_m)
      	        a=a[1..9,.]
            	a=(a-p):/p
	            temp[j,.]=a'
	      }
	bslogit[z,.]=100*47/(colsum(J(47,1,1)-!rowsum(temp)))*mean(temp,1)
	}
end
mata: st_matrix("bslogit", bslogit)
drop _all
svmat bslogit
save "c:\data\restat_12_9\bootstrap_results\post_iv_syruplogit", replace
